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Solar Modulation of Galactic Cosmic Rays, 2 

by 


L. A. Fisk* 

NASA/Goddard Space Flight Center 
Greenbelt, Maryland 


The modulation of galactic cosmic rays in the inter- 
planetary medium can be discussed in terms of a spherically 
symmetric model in which particles undergo convection, dif- 
fusion, and energy changes resulting from the expansion of 
the solar wind. In this model a Fokker-Planck equation de- 
termines, in principle, the particle number density when the 
solar wind speed, the diffusion coefficient, and the inter- 
stellar cosmic ray spectrum are specified (Parker, 1965, 1966; 
Gleeson and Axford, 1967, 1968a, b) . This equation, however, 
is difficult to solve analytically, and, indeed, analytic 
solutions valid at all energies with realistic forms of the 
diffusion coefficient have not been obtained (see the first 
paper in this series, (Fisk and Axford, 1969)). It can be 
solved numerically, and in this Letter we will outline an 
appropriate numerical technique, and present as examples of 
the use of this technique some numerical solutions that pro- 
vide reasonable fits to observed spectra of protons and helium 

♦NAS/NASA Resident Postdoctoral Research Associate 


- 1 - 





- 2 - 


ions with the assumption of quite simple forms for the dif- 
fusion coefficient and unmodulated spectra. 

In the quasi-steady, spherically symmetric model of the 
interplanetary medium discussed by Parker (1965) and by Gleeson 
and Axford (1967), the cosmic ray number density U(r,T), per 
unit interval of kinetic energy T, satisfies a Fokker-Planck 
equat ion : 


? h (r2w > - r-r 
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The corresponding streaming S(r,T) (radial current density), 
per unit interval of kinetic energy, is determined by (Gleeson 
and Axford, 1967): 


s = ™ - 4? - 1 4 ( * TU) 


( 2 ) 


Here K(r,T) is the particle diffusion coefficient, V(r) is 

the solar wind speed, and a(T) = (T+2 T q )/(t + T 0 ) with T c the rest 

energy of a particle. 

Equation (1) is a parabolic partial differential equa- 
tion, and with certain modifications an be solved using the 
numerical techniques which have been developed for dealing 
with simple space-time diffusion equations. Kinetic energy 
in equation (1) is the analogue of time in the simple diffusion 
equation. We will use here the Crank-Nicholson implicit fi- 
nite difference technique (Diaz, 1958), which has the advan- 
tage over other techniques in that in general it determines 
solutions that are unconditionally stable. In this Letter we 
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will outline only the required modifications to equation (1). 

For details on the use of the Crank-Nicholson technique the 
reader is referred to the discussion by Diaz (1958). 

The Crank-Nicholson technique can be used in determin- 
ing the number density in the region 0<r; R, where R is the 
outer boundary of the modulating region, e.g. R could be some 
radial distance where the modulation becomes negligibly small. 

We assume that we know K(r,T), V(r), and also the unmodulated 
spectrum, U(r,T) = U 0 (T) at r = R. Further, we must specify 
a boundary condition at the origin (r = 0) , and also an 
"initial" condition, i„e. , U(r,T) at a given value of T for 
all r , (0- r- R) . In specifying the former condition we could 
require that r 2 S-0 as r-O v i.e., no sources or sinks at the 
origin) since we are concerned here only with galactic cosmic 
rays. However, this condition is difficult to treat numerical- 
ly, and it is considerably simpler merely to note that for re- 
asonable solutions for U, r 2 U-*0 as r-0. We then rewrite equation 
(1) in terms of the variable u(r,T) = r 2 U(r,T): 


!?< Vu> - h 2 l? (rSv) if (oTu) - 1? H?)- 2 !? (?) (3) 

and solve this new equation subject to the boundary conditions 

u(r,T) = 0 at r = 0 , and u(r,T) = R S U 0 (T) at r = R. u(r,T) 
in turn determines U(r,T), but note that we will be unable 
to obtain the number density accurately in the immediate vi- 
cinity of the origin. To specify the "initial" condition 
appropriate for equation (3) (i.e. u(r,T) at a given value 
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of T for all r (0^r*R)), we note that at very large energies 
the effects of modulation can be neglected at all values of 
r. The appropriate initial condition is then u(r,T) = r s U 0 (T) 
at T = T', (0<r- R) , where T' is a sufficiently large energy. 
Clearly, the choice for T' is dependent on the choice for the 
diffusion coefficient, but with realistic values of K , T' ~ 

60 GeV/nucleon should be adequate. With these forms for the 
boundary and "initial" conditions, the Crank-Nicholson tech- 
nique can be used in determining u(r,T) and hence U(r,T), 
beginning with u(r,T) at T = T*, (Ch rcR) , and then calculating 
u(r,T) at lower energies in a step-by-step manner. 

There is an alternative method for specifying the "initial" 
condition appropriate for equation (3). Gleeson and Axford 
(1968b) and Fisk and Axford (1969) have shown that the stream- 
ing S can be neglected in equation (3) when the particles under- 
go relatively little modulation, i.e., when the parameter ft = 
Vr/K^ 1 (here the tilde denotes characteristic value) . The re- 
sulting approximate equation, 

w - I si <aTUfc0 <4) 

can often be easily solved (see Gleeson and Axford, 1968b), 
and its solutions, U a (r,T), used to specify the initial con- 
dition: u(r,T) = r R U a (r,T) at T = T' ' , (0<r- R) where T ~ is 
some energy at which ft^l for all r. In practice, we anticipate 
that at earth ft <1 at energies above a few hundred MeV/nucleon, 
and accordingly equation (4) can be used to determine an ac- 
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curate approximation to the number density at intermediate 
and high energies. Throughout the modulating region (0* r R) , 
at energies above, say, 1 GeV/nucleon. Clearly, the ad- 
vantage in using this "initial" condition over the one de- 
scribed above is that here the computations are begun a much 
lower energy (T " ~ 1 GeV/nucleon vs. T ' ~ 60 GeV/nucleon), 
and hence considerable computational time can be saved in 
determining the modulation at energies T ~ 0.01 - 1.0 GeV/ 
nucleon . 

On assuming forms for U Q , k , and V, the numerical tech- 
nique outlined here can be used to compute the number densi- 
ties and hence the intensities at earth for different species 
of particles, e.g., for protons and helium ions. (Note: in- 
tensity J = vU/4i- , where v is particle speed). These pre- 
dicted intensities can then be compared with observed spectra, 
thereby testing both the assumed forms for U Q , k , and V, and 
also the theory itself. There is, however, considerable lati- 
tude available to us in the choice of the forms for U a and K 
to be used, and hence no definite conclusions will be pos- 
sible without an exhaustive study of the consequences of all 
possible forms. We will present here only one set of pos- 
sible numerical solutions for the intensities of protons and 
helium ions, which we have obtained using simple but realistic 
forms for U Q , < , and V. 

We assume that the unmodulated spectra of both protons 
and helium ions is described by a power law in total energy, 


- 6 - 


-2 65 

U q (T) = A (T+T 0 ) ’ . The constant A is chosen so that 

the unmodulated spectra match the observed spectra at high 
energies where the effects of modulation can be neglected. 

In choosing the form for the diffusion coefficient to be used 
it is important to remember that at a given value of the radial 
distance acceptable diffusion coefficients must be expressible 
in the form particle velocity times a function of rigidity, 
where this function is the same for all species of particles. 
Hence, expressed in the form K = gK-, (p,r) , where g = v/c (c 
is the speed of light) and P is particle rigidity, our choice 
for the diffusion coefficient must be the same for both protons 
and helium ions. We assume here that K = tc o 0Pexp(r/r o ) , where 
K Q is a constant that can be adjusted so that the predicted 
and observed intensities are in the best possible agreement. 

The assumption that ((“SP has been predicted and to some ex- 
tent confirmed (Jokipii, 1966; Gloeckler and Jokipii, 1966; 
O'Gallagher, 1967), and the form K-=exp(r/r 0 ) is a useful one 
since then the dependence of the modulation on radial distance 
is described by a single parameter, a fall-off distance r Q 
which we take to be 1 A.U. in agreement with the findings of 
Gleeson and Axford (1968b) . We take the boundary of the modu- 
lating region to be at R = 25 A.U. (the modulation at this 
point is negligibly small) , and we take V the solar wind speed 
to be constant at 400 km. sec. “1. 

In figure 1 we compare a numerical solution for the inten- 
sity, determined using the above forms for U g ,k, and V, with 
the available observations of proton intensities during solar 
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minimum conditions in 1965. In figure 2 we show the cor- 
responding helium data and numerical solution. k q is taken 
to be k q = 8x10*^ cm.^sec. -1 MV - \ 

e.g., at r = 1 A.U. and T = 6 GeV/nucleon, k = 1.38x10^ 
cm.^sec."^ for protons, and K = 2 . 75 x 10 ^ cm.^sec. - ^ for 
helium ions. As is evident in the figures, the predicted 
intensities agree reasonably well with the observations, 
except at low energies where neither of predicted curves is a 
particularly good fit. 

In figures 1 and 2 we have also plotted the intensities 
determined by equation (4) for the forms for U Q , K , and V 
considered here. As we indicated above, we anticipate that 
at earth equation (4) determines an accurate approximation to 
the number density and hence the intensity at energies above 
a few hundred MeV/nucleon. This equation is known as the 
"force-field" equation since particles behave in this approxi- 
mation as if they were modulated by a heliocentric force field 
(see Gleeson and Axford, 1968b). Indeed, the force-field 
solutions shown in figures 1 and 2 are scarcely distinguishable 
from the corresponding numerical solutions down to energies of 
about 150 MeV/nucleon. In figures 1 and 2 we havfe also shown 
the intensities determined by the simple convection-diffusion 
equation, 

VU = , (5) 

which was the basic equation used in earlier theories that 
neglected the effects of particle energy changes (Parker, 1963). 




- 8 - 


The differences between the simple convection-diffusion 
solutions and the corresponding numerical solutions at 
energies below ~ 400 MeV/nucleon clearly indicate the in- 
adequacies of the earlier simple theory. 

In summary, we have outlined in the Letter a numerical 
technique which can be used to solve the Fokker-Planck 
equation, equation (1), and we have presented some numerical 
solutions that provided reasonable fits to observed spectra of 
protons and helium ions with the assumption of quite simple 
forms of the unmodulated spectra (U c ) and the diffusion co- 
efficient (k). We should not conclude from this, however, 
that we have in fact found the best or even the most likely 
forms for U Q and k since other combinations of these quantities 
will probably yield spectra in better agreement with the 
observations. Rather, the numerical solutions presented 
here should be considered only as examples illustrating the 
use of the numerical technique. 
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FIGURE 1 

A comparison between a numerical solution for the 
proton intensity, and the observations of proton intensi- 
ties during solar minimum conditions in 1965. The inten- 
sities determined by the force-field equation (equation (4)) 
and the simple convection-diffusion equation (equation (5)), 
using the same forms of U Q , K and V as in the numerical solu- 
tion, are also shown. The unmodulated spectrum shown in this 
figure is a plot of the intensity corresponding to the unmodu- 
lated number density; J = vU q /4tt, where U Q is given by a power 
law in total energy. The data points are taken from a paper 
by Gloeckler and Jokipii (1967); the symbol 0 is used to rep- 
resent the observations of Fan, et al. (1966) ; the symbol < 
represents the observations of Balasubrahmanyan , et al. 

(1966a, b) ; the symbol A represents the observations of 
Waddington and Freier (1966); the symbol V represents the 
observations of Ormes and Webber (1966) ; and the symbol > , 
the observations of McDonald (1958). 


FIGURE 2 


A comparison between a numerical solution for the 
helium ion intensity, and the observations of helium ion 
intensities during solar minimum conditions in 1965. The 
intensities determined by the force-field equation (equation 

(4) ) and the simple convection-diffusion equation (equation 

(5) ), using the same forms of U Q , K, and V as in the numerical 
solution, are also shown. The unmodulated spectrum shown in 
this figure is a plot of the intensity corresponding to the 
unmodulated number density; J = vU q /4tt , where U Q is given by 

a power law in total energy. The data points are taken from 
a paper by Gloeckler and Jokipii (1967); the symbol 0 is used 
to represent the observations of Fan, £t al. (1966) ; the 
symbol <i represents the observations of Balasubrahmanyan , 
et al . (1966a, b) ; the symbol A represents the observations of 
Freier and Waddington (1965) ; the symbol v represents the 
observations of Ormes and Webber (1966) ; the symbol 0 re- 
presents the observations of Hofmann and Winckler (1966) ; 
and the symbol o the observations of Anand, et al. (1966). 
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